Multivariate analysis
if (rerun_MCMC) {
HS <- h2life %>%
filter(treat == "HS") %>% rename(Age_HS = NewAge) %>%
as.data.frame()
DR <- h2life %>%
filter(treat == "DR") %>% rename(Age_DR = NewAge) %>%
as.data.frame()
STD <- h2life %>%
filter(treat == "STD") %>% rename(Age_STD = NewAge) %>%
as.data.frame()
h2life_mrg <- full_join(full_join(HS, DR), STD)
prior1 <- list(R = list(V = diag(3) * 1.002, nu = 1.002),
G = list(G1 = list(V = diag(3) * 1.002, nu = 0.002)))
# prior2 <- list(R = list(V = diag(3) / 3, nu = 1.002),
# G = list(G1 = list(V = diag(3) / 3, nu = 0.002)))
#
# prior3 <- list(R = list(V = diag(3) * 1.002, nu = 1),
# G = list(G1 = list(V = diag(3) * 0.5, nu = 0.002)))
# priors <- list(prior1, prior2, prior3)
# prior_names <- c("Tri: V = diag(3) * 1.002, nu = 0.02",
# "Tri: V = diag(3) / 3, nu = 0.02",
# "Tri: V = diag(3) * 0.5, nu = 0.02")
priors <- list(prior1)
prior_names <- c("Tri: V = diag(3) * 1.002, nu = 0.02")
for (ii in 1:length(priors)) {
prior <- priors[[ii]]
fm <- mclapply(1:chains, function(i) {
MCMCglmm(cbind(Age_STD, Age_HS, Age_DR) ~ trait - 1,
random = ~ us(trait):animal,
rcov = ~ idh(trait):units,
data = h2life_mrg,
prior = prior,
pedigree = pedigree,
family = rep("gaussian", 3),
nitt = iter,
burnin = burnin,
thin = thinning)
}, mc.cores = chains)
outfile <- paste0("../../Data/Processed/surv_multivariate_model_prior", ii, ".Rda")
save(fm, file = outfile)
re <- lapply(fm, function(m) m$VCV)
re <- do.call(mcmc.list, re)
re <- as.mcmc(as.matrix(re))
n_eff <- effectiveSize(re)
# STD vs. HS
HS_STD <- re[ , "traitAge_HS:traitAge_STD.animal"] /
sqrt(re[ , "traitAge_STD:traitAge_STD.animal"] *
re[ , "traitAge_HS:traitAge_HS.animal"])
# STD vs. DR
DR_STD <- re[ , "traitAge_DR:traitAge_STD.animal"] /
sqrt(re[ , "traitAge_STD:traitAge_STD.animal"] *
re[ , "traitAge_DR:traitAge_DR.animal"])
# DR vs. HS
HS_DR <- re[ , "traitAge_HS:traitAge_DR.animal"] /
sqrt(re[ , "traitAge_DR:traitAge_DR.animal"] *
re[ , "traitAge_HS:traitAge_HS.animal"])
genet_corr <- bind_rows(genet_corr,
tibble(model = prior_names[[ii]],
HS_STD = median(HS_STD),
DR_STD = median(DR_STD),
HS_DR = median(HS_DR),
n_eff = mean(n_eff)))
}
write_csv(genet_corr, path = "../../Data/Processed/Lifespan_Genetic_Correlations.csv")
}
Analyze model
load("../../Data/Processed/surv_multivariate_model_prior1.Rda")
fe <- lapply(fm, function(m) m$Sol)
fe <- do.call(mcmc.list, fe)
plot(fe[, 1, drop = FALSE], ask = FALSE)

plot(fe[, 2, drop = FALSE], ask = FALSE)

plot(fe[, 3, drop = FALSE], ask = FALSE)

# Extract random effects, convert to mcmc.list
re <- lapply(fm, function(m) m$VCV)
re <- do.call(mcmc.list, re)
plot(re[, 1, drop = FALSE], ask = FALSE)

plot(re[, 2, drop = FALSE], ask = FALSE)

plot(re[, 3, drop = FALSE], ask = FALSE)

autocorr.diag(re)
## traitAge_STD:traitAge_STD.animal traitAge_HS:traitAge_STD.animal
## Lag 0 1.000000000 1.0000000000
## Lag 500 0.068179841 0.0736391590
## Lag 2500 -0.006882237 -0.0038467280
## Lag 5000 0.007044146 -0.0061371156
## Lag 25000 0.004727290 0.0009293054
## traitAge_DR:traitAge_STD.animal traitAge_STD:traitAge_HS.animal
## Lag 0 1.000000000 1.0000000000
## Lag 500 0.063500507 0.0736391590
## Lag 2500 -0.004802190 -0.0038467280
## Lag 5000 -0.001234391 -0.0061371156
## Lag 25000 -0.001749820 0.0009293054
## traitAge_HS:traitAge_HS.animal traitAge_DR:traitAge_HS.animal
## Lag 0 1.000000000 1.0000000000
## Lag 500 0.083782065 0.0888369970
## Lag 2500 -0.017718219 0.0045442414
## Lag 5000 -0.005514260 -0.0107881960
## Lag 25000 -0.001397025 -0.0004649028
## traitAge_STD:traitAge_DR.animal traitAge_HS:traitAge_DR.animal
## Lag 0 1.000000000 1.0000000000
## Lag 500 0.063500507 0.0888369970
## Lag 2500 -0.004802190 0.0045442414
## Lag 5000 -0.001234391 -0.0107881960
## Lag 25000 -0.001749820 -0.0004649028
## traitAge_DR:traitAge_DR.animal traitAge_STD.units
## Lag 0 1.000000000 1.000000000
## Lag 500 0.085607455 0.059879186
## Lag 2500 0.002057077 -0.007839362
## Lag 5000 0.004328985 -0.001042706
## Lag 25000 0.006550331 0.001383093
## traitAge_HS.units traitAge_DR.units
## Lag 0 1.0000000000 1.000000000
## Lag 500 0.0665168077 0.062501759
## Lag 2500 -0.0165415030 0.008408074
## Lag 5000 -0.0014701485 -0.005139524
## Lag 25000 -0.0007576496 0.008627476
effectiveSize(re)
## traitAge_STD:traitAge_STD.animal traitAge_HS:traitAge_STD.animal
## 21514.52 20212.93
## traitAge_DR:traitAge_STD.animal traitAge_STD:traitAge_HS.animal
## 22013.61 20212.93
## traitAge_HS:traitAge_HS.animal traitAge_DR:traitAge_HS.animal
## 19734.78 19002.38
## traitAge_STD:traitAge_DR.animal traitAge_HS:traitAge_DR.animal
## 22013.61 19002.38
## traitAge_DR:traitAge_DR.animal traitAge_STD.units
## 20088.25 21512.98
## traitAge_HS.units traitAge_DR.units
## 20320.53 20642.51
# Concatenate samples
re <- as.mcmc(as.matrix(re))
head(re)
## Markov Chain Monte Carlo (MCMC) output:
## Start = 1
## End = 7
## Thinning interval = 1
## traitAge_STD:traitAge_STD.animal traitAge_HS:traitAge_STD.animal
## [1,] 0.4632983 0.16953957
## [2,] 0.3950603 0.13623892
## [3,] 0.3263196 0.01576263
## [4,] 0.5508081 0.12743270
## [5,] 0.5017048 0.12533453
## [6,] 0.5403548 0.25806069
## [7,] 0.3794987 0.10694896
## traitAge_DR:traitAge_STD.animal traitAge_STD:traitAge_HS.animal
## [1,] 0.1009902 0.16953957
## [2,] 0.1959820 0.13623892
## [3,] 0.1045252 0.01576263
## [4,] 0.1523802 0.12743270
## [5,] 0.1975044 0.12533453
## [6,] 0.1772361 0.25806069
## [7,] 0.1167529 0.10694896
## traitAge_HS:traitAge_HS.animal traitAge_DR:traitAge_HS.animal
## [1,] 0.2790736 0.11957617
## [2,] 0.3379790 0.10244792
## [3,] 0.2690235 0.07564472
## [4,] 0.3205873 0.06663020
## [5,] 0.3030006 0.07402833
## [6,] 0.3937808 0.13460634
## [7,] 0.3468938 0.12529748
## traitAge_STD:traitAge_DR.animal traitAge_HS:traitAge_DR.animal
## [1,] 0.1009902 0.11957617
## [2,] 0.1959820 0.10244792
## [3,] 0.1045252 0.07564472
## [4,] 0.1523802 0.06663020
## [5,] 0.1975044 0.07402833
## [6,] 0.1772361 0.13460634
## [7,] 0.1167529 0.12529748
## traitAge_DR:traitAge_DR.animal traitAge_STD.units traitAge_HS.units
## [1,] 0.2258080 0.4949568 0.5178935
## [2,] 0.2310788 0.5079044 0.4634037
## [3,] 0.1762062 0.5622988 0.5508867
## [4,] 0.1506864 0.4227779 0.4810127
## [5,] 0.2671842 0.4244858 0.4844486
## [6,] 0.1886153 0.4402332 0.4419436
## [7,] 0.2341845 0.5260786 0.4923344
## traitAge_DR.units
## [1,] 0.5405749
## [2,] 0.5475709
## [3,] 0.5945725
## [4,] 0.6132244
## [5,] 0.5262783
## [6,] 0.5440586
## [7,] 0.5453158
# STD vs. HS
plot(re[ , "traitAge_HS:traitAge_STD.animal"])

plot(re[ , "traitAge_STD:traitAge_STD.animal"])

HS_STD <- re[ , "traitAge_HS:traitAge_STD.animal"] /
sqrt(re[ , "traitAge_STD:traitAge_STD.animal"] *
re[ , "traitAge_HS:traitAge_HS.animal"])
plot(HS_STD)

median(HS_STD)
## [1] 0.2999119
HPDinterval(HS_STD)
## lower upper
## var1 0.02156968 0.5526564
## attr(,"Probability")
## [1] 0.95
# STD vs. DR
DR_STD <- re[ , "traitAge_DR:traitAge_STD.animal"] /
sqrt(re[ , "traitAge_STD:traitAge_STD.animal"] *
re[ , "traitAge_DR:traitAge_DR.animal"])
plot(DR_STD)

median(DR_STD)
## [1] 0.4957338
HPDinterval(DR_STD)
## lower upper
## var1 0.2584886 0.7007605
## attr(,"Probability")
## [1] 0.95
# DR vs. HS
HS_DR <- re[ , "traitAge_HS:traitAge_DR.animal"] /
sqrt(re[ , "traitAge_DR:traitAge_DR.animal"] *
re[ , "traitAge_HS:traitAge_HS.animal"])
plot(HS_DR)

median(HS_DR)
## [1] 0.3501346
HPDinterval(HS_DR)
## lower upper
## var1 0.07818288 0.6078474
## attr(,"Probability")
## [1] 0.95
Plot for paper
M <- data_frame(`HS vs. STD` = as.numeric(HS_STD),
`DR vs. STD` = as.numeric(DR_STD),
`HS vs. DR` = as.numeric(HS_DR))
M <- as_tibble( M ) %>%
select(`HS vs. STD`, `DR vs. STD`, `HS vs. DR`) %>%
rename(`HS vs. C` = `HS vs. STD`) %>%
rename(`DR vs. C` = `DR vs. STD`) %>%
rename(`HS vs. DR` = `HS vs. DR`)
colMeans(M)
## HS vs. C DR vs. C HS vs. DR
## 0.2939628 0.4873072 0.3441595
M %>% gather(Comparison, value) %>%
ggplot(aes(value, color = Comparison)) +
geom_line(stat = "density", size = 1.0) +
labs(x = "Genetic Correlation", y = "Density") +
theme(legend.position = c(0.10, 0.85),
text = element_text(size = 10),
legend.text = element_text(size = 10)) +
scale_x_continuous(limits = c(-1, 1)) +
my_theme

ggsave(file = "../../Figures/Lifespan_Genetic_Correlation_Plot.pdf",
width = 4, height = 4)
Lifespan_correlation <- last_plot()
save(Lifespan_correlation, file = "../../Figures/Lifespan_correlation.Rda")
Parameter expanded prior multivariate analysis
if (rerun_MCMC) {
iter <- 1e5
burnin <- 2e3
thinning <- 500
chains <- 12
HS <- h2life %>%
filter(treat == "HS") %>% rename(Age_HS = NewAge) %>%
as.data.frame()
DR <- h2life %>%
filter(treat == "DR") %>% rename(Age_DR = NewAge) %>%
as.data.frame()
STD <- h2life %>%
filter(treat == "STD") %>% rename(Age_STD = NewAge) %>%
as.data.frame()
h2life_mrg <- full_join(full_join(HS, DR), STD)
prior1 <- list(R = list(V = diag(3), nu = 1.002, fix = 3),
G = list(G1 = list(V = diag(3),
nu = 1,
alpha.mu = rep(0, 3),
alpha.V = diag(3) * 25 ^ 2)))
priors <- list(prior1)
prior_names <- c("Tri: Parameter expanded")
for (ii in 1:length(priors)) {
prior <- priors[[ii]]
fm <- mclapply(1:chains, function(i) {
MCMCglmm(cbind(Age_STD, Age_HS, Age_DR) ~ trait - 1,
random = ~ us(trait):animal,
rcov = ~ idh(trait):units,
data = h2life_mrg,
prior = prior,
pedigree = pedigree,
family = rep("gaussian", 3),
nitt = iter,
burnin = burnin,
thin = thinning)
}, mc.cores = chains)
outfile <- paste0("../../Data/Processed/surv_multivariate_model_prior_expanded", ii, ".Rda")
save(fm, file = outfile)
re <- lapply(fm, function(m) m$VCV)
re <- do.call(mcmc.list, re)
re <- as.mcmc(as.matrix(re))
n_eff <- effectiveSize(re)
# STD vs. HS
HS_STD <- re[ , "traitAge_HS:traitAge_STD.animal"] /
sqrt(re[ , "traitAge_STD:traitAge_STD.animal"] *
re[ , "traitAge_HS:traitAge_HS.animal"])
# STD vs. DR
DR_STD <- re[ , "traitAge_DR:traitAge_STD.animal"] /
sqrt(re[ , "traitAge_STD:traitAge_STD.animal"] *
re[ , "traitAge_DR:traitAge_DR.animal"])
# DR vs. HS
HS_DR <- re[ , "traitAge_HS:traitAge_DR.animal"] /
sqrt(re[ , "traitAge_DR:traitAge_DR.animal"] *
re[ , "traitAge_HS:traitAge_HS.animal"])
genet_corr <- bind_rows(genet_corr,
tibble(model = prior_names[[ii]],
HS_STD = median(HS_STD),
DR_STD = median(DR_STD),
HS_DR = median(HS_DR),
n_eff = mean(n_eff)))
}
write_csv(
genet_corr,
path = "../../Data/Processed/Lifespan_Genetic_Correlations_Parameter_Expanded.csv")
}
Analyze parameter expanded model
load("../../Data/Processed/surv_multivariate_model_prior_expanded1.Rda")
fe <- lapply(fm, function(m) m$Sol)
fe <- do.call(mcmc.list, fe)
plot(fe[, 1, drop = FALSE], ask = FALSE)

plot(fe[, 2, drop = FALSE], ask = FALSE)

plot(fe[, 3, drop = FALSE], ask = FALSE)

# Extract random effects, convert to mcmc.list
re <- lapply(fm, function(m) m$VCV)
re <- do.call(mcmc.list, re)
plot(re[, 1, drop = FALSE], ask = FALSE)

plot(re[, 2, drop = FALSE], ask = FALSE)

plot(re[, 3, drop = FALSE], ask = FALSE)

autocorr.diag(re)
## traitAge_STD:traitAge_STD.animal traitAge_HS:traitAge_STD.animal
## Lag 0 1.000000000 1.000000000
## Lag 500 -0.011880172 0.084936923
## Lag 2500 -0.012336107 -0.044213068
## Lag 5000 -0.008184613 0.022182921
## Lag 25000 0.009565219 0.005296313
## traitAge_DR:traitAge_STD.animal traitAge_STD:traitAge_HS.animal
## Lag 0 1.000000000 1.000000000
## Lag 500 0.077587962 0.084936923
## Lag 2500 -0.010870186 -0.044213068
## Lag 5000 0.030051023 0.022182921
## Lag 25000 -0.007022395 0.005296313
## traitAge_HS:traitAge_HS.animal traitAge_DR:traitAge_HS.animal
## Lag 0 1.0000000000 1.000000000
## Lag 500 0.0003537141 0.014058983
## Lag 2500 0.0136531239 -0.023536549
## Lag 5000 -0.0173873268 -0.017871792
## Lag 25000 -0.0131964670 0.009574731
## traitAge_STD:traitAge_DR.animal traitAge_HS:traitAge_DR.animal
## Lag 0 1.000000000 1.000000000
## Lag 500 0.077587962 0.014058983
## Lag 2500 -0.010870186 -0.023536549
## Lag 5000 0.030051023 -0.017871792
## Lag 25000 -0.007022395 0.009574731
## traitAge_DR:traitAge_DR.animal traitAge_STD.units
## Lag 0 1.000000000 1.000000000
## Lag 500 0.003887223 -0.002186594
## Lag 2500 -0.005213611 -0.024391116
## Lag 5000 -0.033258672 -0.009846700
## Lag 25000 0.001597150 0.006159378
## traitAge_HS.units traitAge_DR.units
## Lag 0 1.000000000 NaN
## Lag 500 0.023515392 NaN
## Lag 2500 -0.001714462 NaN
## Lag 5000 -0.036410907 NaN
## Lag 25000 -0.008078434 NaN
effectiveSize(re)
## traitAge_STD:traitAge_STD.animal traitAge_HS:traitAge_STD.animal
## 2339.831 2209.375
## traitAge_DR:traitAge_STD.animal traitAge_STD:traitAge_HS.animal
## 2331.718 2209.375
## traitAge_HS:traitAge_HS.animal traitAge_DR:traitAge_HS.animal
## 2737.938 2323.729
## traitAge_STD:traitAge_DR.animal traitAge_HS:traitAge_DR.animal
## 2331.718 2323.729
## traitAge_DR:traitAge_DR.animal traitAge_STD.units
## 2475.675 2512.021
## traitAge_HS.units traitAge_DR.units
## 2437.652 0.000
# Concatenate samples
re <- as.mcmc(as.matrix(re))
head(re)
## Markov Chain Monte Carlo (MCMC) output:
## Start = 1
## End = 7
## Thinning interval = 1
## traitAge_STD:traitAge_STD.animal traitAge_HS:traitAge_STD.animal
## [1,] 0.3177375 0.0001236148
## [2,] 0.5171791 0.0012672359
## [3,] 0.3841529 0.0197225670
## [4,] 0.3423261 0.0346445797
## [5,] 0.3976341 0.0676707225
## [6,] 0.5314955 0.1152669886
## [7,] 0.3384343 0.1062301494
## traitAge_DR:traitAge_STD.animal traitAge_STD:traitAge_HS.animal
## [1,] 0.08665067 0.0001236148
## [2,] 0.09078737 0.0012672359
## [3,] 0.07710740 0.0197225670
## [4,] 0.04821795 0.0346445797
## [5,] 0.06409478 0.0676707225
## [6,] 0.06619153 0.1152669886
## [7,] 0.05919508 0.1062301494
## traitAge_HS:traitAge_HS.animal traitAge_DR:traitAge_HS.animal
## [1,] 0.2019285 0.03867934
## [2,] 0.3023297 0.06524771
## [3,] 0.2444371 0.07475899
## [4,] 0.2704081 0.06776188
## [5,] 0.2260605 0.06126987
## [6,] 0.2762569 0.08347054
## [7,] 0.2542423 0.07237181
## traitAge_STD:traitAge_DR.animal traitAge_HS:traitAge_DR.animal
## [1,] 0.08665067 0.03867934
## [2,] 0.09078737 0.06524771
## [3,] 0.07710740 0.07475899
## [4,] 0.04821795 0.06776188
## [5,] 0.06409478 0.06126987
## [6,] 0.06619153 0.08347054
## [7,] 0.05919508 0.07237181
## traitAge_DR:traitAge_DR.animal traitAge_STD.units traitAge_HS.units
## [1,] 0.07656928 0.5323402 0.5340417
## [2,] 0.07191132 0.4571359 0.4873783
## [3,] 0.08480175 0.4993428 0.5207448
## [4,] 0.05734316 0.5560712 0.4897561
## [5,] 0.07447087 0.5157063 0.5407287
## [6,] 0.09450949 0.4034260 0.5127812
## [7,] 0.10054775 0.5484401 0.5123482
## traitAge_DR.units
## [1,] 1
## [2,] 1
## [3,] 1
## [4,] 1
## [5,] 1
## [6,] 1
## [7,] 1
# STD vs. HS
plot(re[ , "traitAge_HS:traitAge_STD.animal"])

plot(re[ , "traitAge_STD:traitAge_STD.animal"])

HS_STD <- re[ , "traitAge_HS:traitAge_STD.animal"] /
sqrt(re[ , "traitAge_STD:traitAge_STD.animal"] *
re[ , "traitAge_HS:traitAge_HS.animal"])
plot(HS_STD)

median(HS_STD)
## [1] 0.2311976
HPDinterval(HS_STD)
## lower upper
## var1 -0.03035282 0.5091342
## attr(,"Probability")
## [1] 0.9498299
# STD vs. DR
DR_STD <- re[ , "traitAge_DR:traitAge_STD.animal"] /
sqrt(re[ , "traitAge_STD:traitAge_STD.animal"] *
re[ , "traitAge_DR:traitAge_DR.animal"])
plot(DR_STD)

median(DR_STD)
## [1] 0.4468741
HPDinterval(DR_STD)
## lower upper
## var1 0.2074835 0.6598369
## attr(,"Probability")
## [1] 0.9498299
# DR vs. HS
HS_DR <- re[ , "traitAge_HS:traitAge_DR.animal"] /
sqrt(re[ , "traitAge_DR:traitAge_DR.animal"] *
re[ , "traitAge_HS:traitAge_HS.animal"])
plot(HS_DR)

median(HS_DR)
## [1] 0.3498687
HPDinterval(HS_DR)
## lower upper
## var1 0.112513 0.6110439
## attr(,"Probability")
## [1] 0.9498299
Plot parameter expanded model
M <- data_frame(`HS vs. STD` = as.numeric(HS_STD),
`DR vs. STD` = as.numeric(DR_STD),
`HS vs. DR` = as.numeric(HS_DR))
M <- as_tibble( M ) %>%
select(`HS vs. STD`, `DR vs. STD`, `HS vs. DR`) %>%
rename(`HS vs. C` = `HS vs. STD`) %>%
rename(`DR vs. C` = `DR vs. STD`) %>%
rename(`HS vs. DR` = `HS vs. DR`)
colMeans(M)
## HS vs. C DR vs. C HS vs. DR
## 0.2304897 0.4403069 0.3474292
M %>% gather(Comparison, value) %>%
ggplot(aes(value, color = Comparison)) +
geom_line(stat = "density", size = 1.0) +
labs(x = "Genetic Correlation", y = "Density") +
theme(legend.position = c(0.10, 0.85),
text = element_text(size = 10),
legend.text = element_text(size = 10)) +
scale_x_continuous(limits = c(-1, 1)) +
my_theme

# ggsave(file = "../../Figures/Lifespan_Genetic_Correlation_Plot.pdf",
# width = 4, height = 4)
# Lifespan_correlation <- last_plot()
# save(Lifespan_correlation, file = "../../Figures/Lifespan_correlation.Rda")